# Script to produce Table 1, Figure 1
# Gould et al. 2023 EHP

# 0. Libraries -------

library(tidyverse)
library(haven)
library(cowplot)
library(ggstatsplot)
library(broom)

# Specify how we want summary tables to look (from arsenal package)
mycontrols  <- tableby.control(test=TRUE, total=FALSE,
                               numeric.stats=c("meansd", "medianq1q3"),
                               cat.stats=c("countpct"),
                               stats.labels=list(countpct='N (%)', meansd='Mean (SD)', medianq1q3="Median (IQR)", range = "Range"),
                               digits = 2L,
                               digits.pct = 0L,
                               digits.p = 3L,
                               pfootnote = T)

# 1. Data ---------

# Main data ------
full_df <- read_rds("~/Downloads/gould_ehp_ecuador_cf_df.rds")

# Adult literacy ------

literate_alri <- ggscatterstats(
  data = full_df,
  x = Literate,
  y = ALRI5_5_plus_100kpop,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Adult Literacy (%)", # label for x axis
  ylab = "under-5 LRI mortality\n(per 100,000 under-5 population)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Adult Literacy", # title text for the plot

  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2,
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(),
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    scale_y_log10(breaks = c(0.1, 1, 10, 50, 100, 500, 1000)),
    geom_point(data = full_df,
                    aes(x=Literate, y = ALRI5_5_plus_100kpop,
                        group=period, color=period), alpha=0.2)
    )
)

literate_cf <- ggscatterstats(
  data = full_df,
  x = Literate,
  y = cf,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Adult Literacy (%)", # label for x axis
  ylab = "Clean fuel (%)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Adult Literacy", # title text for the plot

  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2,
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(),
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    ggplot2::scale_y_continuous(labels=scales::percent_format(),
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    geom_point(data = full_df,
               aes(x=Literate, y=cf,
                   group=period, color=period), alpha=0.4)
  )
)

literate_alri_mod <- glm(ALRI5_5_plus ~ offset(log(under5population)) + 
                           Literate + period + canton_id_universal, 
                         full_df,  family=quasipoisson())

literate_cf_mod <- lm(cf ~ Literate + period + canton_id_universal, full_df)

broom::tidy(literate_alri_mod)[2,]
broom::tidy(literate_cf_mod)[2,]

# ggsave("~/Dropbox/Ecuador National Health and LPG Analysis/Final_Repository/Figures/Confounding_Correlations_June1/literate.png",
#        cowplot::plot_grid(
#          literate_alri,
#          literate_cf,
#          align="hv",
#          nrow=2
#        ),
#        width = 220,
#        height = 200,
#        units = "mm", 
#        dpi = 125)



# Adult women literacy ------

women_literate_alri <- ggscatterstats(
  data = full_df,
  x = Literate_women,
  y = ALRI5_5_plus_100kpop,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Adult Women Literacy (%)", # label for x axis
  ylab = "under-5 LRI mortality\n(per 100,000 under-5 population)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Adult Women Literacy", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    scale_y_log10(breaks = c(0.1, 1, 10, 50, 100, 500, 1000)),
    geom_point(data = full_df, 
               aes(x=Literate_women, y = ALRI5_5_plus_100kpop, 
                   group=period, color=period), alpha=0.5)
  )
)

women_literate_cf <- ggscatterstats(
  data = full_df,
  x = Literate_women,
  y = cf,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Adult Women Literacy (%)", # label for x axis
  ylab = "Clean fuel (%)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Adult Women Literacy", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    ggplot2::scale_y_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    geom_point(data = full_df, 
               aes(x=Literate_women, y=cf, 
                   group=period, color=period), alpha=0.4)
  )
)

women_literate_alri
women_literate_cf

literate_women_alri_mod <- glm(ALRI5_5_plus ~ offset(log(under5population)) + 
                                                 Literate_women + period + canton_id_universal, 
                                               full_df,  family=quasipoisson())
literate_women_cf_mod <- lm(cf ~ Literate_women + period + canton_id_universal, full_df)

broom::tidy(literate_women_alri_mod)[2,]
broom::tidy(literate_women_cf_mod)[2,]

# Rural ------

rural_alri <- ggscatterstats(
  data = full_df,
  x = percent_rural_corrected_ratio,
  y = ALRI5_5_plus_100kpop,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Rural (%)", # label for x axis
  ylab = "under-5 LRI mortality\n(per 100,000 under-5 population)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Rural", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    scale_y_log10(breaks = c(0.1, 1, 10, 50, 100, 500, 1000)),
    geom_point(data = full_df, 
               aes(x=percent_rural_corrected_ratio, y=ALRI5_5_plus_100kpop, 
                   group=period, color=period), alpha=0.5)
  )
)

rural_cf <- ggscatterstats(
  data = full_df,
  x = percent_rural_corrected_ratio,
  y = cf,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Rural (%)", # label for x axis
  ylab = "Clean fuel (%)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Rural", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    ggplot2::scale_y_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    geom_point(data = full_df, 
               aes(x=percent_rural_corrected_ratio, y=cf, 
                   group=period, color=period), alpha=0.4)
  )
)

rural_alri
rural_cf

rural_alri_mod <- glm(ALRI5_5_plus ~ offset(log(under5population)) + 
                                 percent_rural_corrected_ratio + period + canton_id_universal, 
                               full_df,  family=quasipoisson())

rural_cf_mod <- lm(cf ~ percent_rural_corrected_ratio + period + canton_id_universal, full_df)

broom::tidy(rural_alri_mod)[2,]
broom::tidy(rural_cf_mod)[2,]



# In school ------

school_alri <- ggscatterstats(
  data = full_df,
  x = InSchool,
  y = ALRI5_5_plus_100kpop,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Children under 18 in school (%)", # label for x axis
  ylab = "under-5 LRI mortality\n(per 100,000 under-5 population)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "In School", # title text for the plot

  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2,
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(),
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    scale_y_log10(breaks = c(0.1, 1, 10, 50, 100, 500, 1000)),
    geom_point(data = full_df,
               aes(x=InSchool, y = ALRI5_5_plus_100kpop,
                   group=period, color=period), alpha=0.5)
  )
)

school_cf <- ggscatterstats(
  data = full_df,
  x = InSchool,
  y = cf,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Children under 18 in school (%)", # label for x axis
  ylab = "Clean fuel (%)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "In School", # title text for the plot

  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2,
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(),
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    ggplot2::scale_y_continuous(labels=scales::percent_format(),
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    geom_point(data = full_df,
               aes(x=InSchool, y=cf,
                   group=period, color=period), alpha=0.4)
  )
)

school_alri
school_cf

# Girls In school ------

girls_school_alri <- ggscatterstats(
  data = full_df,
  x = InSchool_girls,
  y = ALRI5_5_plus_100kpop,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Girls under 18 in school (%)", # label for x axis
  ylab = "under-5 LRI mortality\n(per 100,000 under-5 population)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Girls In School", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    scale_y_log10(breaks = c(0.1, 1, 10, 50, 100, 500, 1000)),
    geom_point(data = full_df, 
               aes(x=InSchool_girls, y = ALRI5_5_plus_100kpop, 
                   group=period, color=period), alpha=0.5)
  )
)

girls_school_cf <- ggscatterstats(
  data = full_df,
  x = InSchool_girls,
  y = cf,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Girls under 18 in school (%)", # label for x axis
  ylab = "Clean fuel (%)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Girls In School", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    ggplot2::scale_y_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    geom_point(data = full_df, 
               aes(x=InSchool_girls, y=cf, 
                   group=period, color=period), alpha=0.4)
  )
)

girls_school_alri
girls_school_cf

girls_school_alri_mod <- glm(ALRI5_5_plus ~ offset(log(under5population)) + 
                               InSchool_girls + period + canton_id_universal, 
                      full_df,  family=quasipoisson())

girls_school_cf_mod <- lm(cf ~ InSchool_girls + period + canton_id_universal, full_df)

broom::tidy(girls_school_alri_mod)[2,]
broom::tidy(girls_school_cf_mod)[2,]


# Basura: service ------

basura_service_alri <- ggscatterstats(
  data = full_df,
  x = elimina_basura_service,
  y = ALRI5_5_plus_100kpop,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Trash removal: service (%)", # label for x axis
  ylab = "under-5 LRI mortality\n(per 100,000 under-5 population)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Trash removal: service", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    scale_y_log10(breaks = c(0.1, 1, 10, 50, 100, 500, 1000)),
    geom_point(data = full_df, 
               aes(x=elimina_basura_service, y = ALRI5_5_plus_100kpop, 
                   group=period, color=period), alpha=0.5)
  )
)

basura_service_cf <- ggscatterstats(
  data = full_df,
  x = elimina_basura_service,
  y = cf,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Trash removal: service (%)", # label for x axis
  ylab = "Clean fuel (%)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Trash removal: service", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    ggplot2::scale_y_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    geom_point(data = full_df, 
               aes(x=elimina_basura_service, y=cf, 
                   group=period, color=period), alpha=0.4)
  )
)

elimina_basura_service_alri_mod <- glm(ALRI5_5_plus ~ offset(log(under5population)) + 
                               elimina_basura_service + period + canton_id_universal, 
                             full_df,  family=quasipoisson())

elimina_basura_service_cf_mod <- lm(cf ~ elimina_basura_service + period + canton_id_universal, full_df)

broom::tidy(elimina_basura_service_alri_mod)[2,]
broom::tidy(elimina_basura_service_cf_mod)[2,]


# Ducha ------

servicio_ducha_alri <- ggscatterstats(
  data = full_df,
  x = servicio_ducha_excluso,
  y = ALRI5_5_plus_100kpop,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Has exclusive shower access (%)", # label for x axis
  ylab = "under-5 LRI mortality\n(per 100,000 under-5 population)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Has exclusive shower access", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    scale_y_log10(breaks = c(0.1, 1, 10, 50, 100, 500, 1000)),
    geom_point(data = full_df, 
               aes(x=servicio_ducha_excluso, y = ALRI5_5_plus_100kpop, 
                   group=period, color=period), alpha=0.5)
  )
)

servicio_ducha_cf <- ggscatterstats(
  data = full_df,
  x = servicio_ducha_excluso,
  y = cf,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Has exclusive shower access (%)", # label for x axis
  ylab = "Clean fuel (%)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Has exclusive shower access", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    ggplot2::scale_y_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    geom_point(data = full_df, 
               aes(x=servicio_ducha_excluso, y=cf, 
                   group=period, color=period), alpha=0.4)
  )
)

servicio_ducha_alri
servicio_ducha_cf

servicio_ducha_alri_mod <- glm(ALRI5_5_plus ~ offset(log(under5population)) + 
                                 servicio_ducha_excluso + period + canton_id_universal, 
                                       full_df,  family=quasipoisson())

servicio_ducha_cf_mod <- lm(cf ~ servicio_ducha_excluso + period + canton_id_universal, full_df)

broom::tidy(servicio_ducha_alri_mod)[2,]
broom::tidy(servicio_ducha_cf_mod)[2,]

# higenico: red publica ------

servicio_higenico_alri <- ggscatterstats(
  data = full_df,
  x = servicio_higenico_redpublica_use,
  y = ALRI5_5_plus_100kpop,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Has modern toilet with municipal system waste removal (%)", # label for x axis
  ylab = "under-5 LRI mortality\n(per 100,000 under-5 population)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Has modern toilet with municipal system waste removal", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    scale_y_log10(breaks = c(0.1, 1, 10, 50, 100, 500, 1000)),
    geom_point(data = full_df, 
               aes(x=servicio_higenico_redpublica_use, y = ALRI5_5_plus_100kpop, 
                   group=period, color=period), alpha=0.5)
  )
)

servicio_higenico_cf <- ggscatterstats(
  data = full_df,
  x = servicio_higenico_redpublica_use,
  y = cf,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Has modern toilet with municipal system waste removal (%)", # label for x axis
  ylab = "Clean fuel (%)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Has modern toilet with municipal system waste removal", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    ggplot2::scale_y_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    geom_point(data = full_df, 
               aes(x=servicio_higenico_redpublica_use, y=cf, 
                   group=period, color=period), alpha=0.4)
  )
)

servicio_higenico_alri
servicio_higenico_cf

servicio_higenico_alri_mod <- glm(ALRI5_5_plus ~ offset(log(under5population)) + 
                                 servicio_higenico_redpublica_use + period + canton_id_universal, 
                               full_df,  family=quasipoisson())

servicio_higenico_cf_mod <- lm(cf ~ servicio_higenico_redpublica_use + period + canton_id_universal, full_df)

broom::tidy(servicio_higenico_alri_mod)[2,]
broom::tidy(servicio_higenico_cf_mod)[2,]

# higenico: red publica, pozo ------

elimina_servidas_red_pozo_alri <- ggscatterstats(
  data = full_df,
  x = elimina_servidas_red_pozo,
  y = ALRI5_5_plus_100kpop,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Has modern toilet with municipal system waste removal, cesspool, or septic tank (%)", # label for x axis
  ylab = "under-5 LRI mortality\n(per 100,000 under-5 population)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Modern toilet and waste removal", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    scale_y_log10(breaks = c(0.1, 1, 10, 50, 100, 500, 1000)),
    geom_point(data = full_df, 
               aes(x=elimina_servidas_red_pozo, y = ALRI5_5_plus_100kpop, 
                   group=period, color=period), alpha=0.5)
  )
)

elimina_servidas_red_pozo_cf <- ggscatterstats(
  data = full_df,
  x = elimina_servidas_red_pozo,
  y = cf,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Has modern toilet with municipal system waste removal, cesspool, or septic tank (%)", # label for x axis
  ylab = "Clean fuel (%)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Modern toilet and waste removal", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    ggplot2::scale_y_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    geom_point(data = full_df, 
               aes(x=elimina_servidas_red_pozo, y=cf, 
                   group=period, color=period), alpha=0.4)
  )
)

elimina_servidas_red_pozo_alri
elimina_servidas_red_pozo_cf

elimina_servidas_red_pozo_alri_mod <- glm(ALRI5_5_plus ~ offset(log(under5population)) + 
                                    elimina_servidas_red_pozo + period + canton_id_universal, 
                                  full_df,  family=quasipoisson())

elimina_servidas_red_pozo_cf_mod <- lm(cf ~ elimina_servidas_red_pozo + period + canton_id_universal, full_df)

broom::tidy(elimina_servidas_red_pozo_alri_mod)[2,]
broom::tidy(elimina_servidas_red_pozo_cf_mod)[2,]

# higenico: red publica, pozo, letrina ------

elimina_servidas_red_pozo_letrina_alri <- ggscatterstats(
  data = full_df,
  x = elimina_servidas_red_pozo_letrina,
  y = ALRI5_5_plus_100kpop,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Has modern toilet with municipal system waste removal, cesspool, septic tank, or latrine (%)", # label for x axis
  ylab = "under-5 LRI mortality\n(per 100,000 under-5 population)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Modern toilet and waste removal", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    scale_y_log10(breaks = c(0.1, 1, 10, 50, 100, 500, 1000)),
    geom_point(data = full_df, 
               aes(x=elimina_servidas_red_pozo_letrina, y = ALRI5_5_plus_100kpop, 
                   group=period, color=period), alpha=0.5)
  )
)

elimina_servidas_red_pozo_letrina_cf <- ggscatterstats(
  data = full_df,
  x = elimina_servidas_red_pozo_letrina,
  y = cf,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Has modern toilet with municipal system waste removal, cesspool, septic tank, or latrine (%)", # label for x axis
  ylab = "Clean fuel (%)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Modern toilet and waste removal", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    ggplot2::scale_y_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    geom_point(data = full_df, 
               aes(x=elimina_servidas_red_pozo_letrina, y=cf, 
                   group=period, color=period), alpha=0.4)
  )
)

elimina_servidas_red_pozo_letrina_alri
elimina_servidas_red_pozo_letrina_cf

elimina_servidas_red_pozo_letrina_alri_mod <- glm(ALRI5_5_plus ~ offset(log(under5population)) + 
                                    elimina_servidas_red_pozo_letrina + period + canton_id_universal, 
                                  full_df,  family=quasipoisson())

elimina_servidas_red_pozo_letrina_cf_mod <- lm(cf ~ elimina_servidas_red_pozo_letrina + period + canton_id_universal, full_df)

broom::tidy(elimina_servidas_red_pozo_letrina_alri_mod)[2,]
broom::tidy(elimina_servidas_red_pozo_letrina_cf_mod)[2,]

# obtiene agua: tuberia ------

agua_tuberia_alri <- ggscatterstats(
  data = full_df,
  x = obtiene_agua_redpublica_tuberia_adentro_use,
  y = ALRI5_5_plus_100kpop,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Gets water via public system through pipes indoors (%)", # label for x axis
  ylab = "under-5 LRI mortality\n(per 100,000 under-5 population)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Piped water", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    scale_y_log10(breaks = c(0.1, 1, 10, 50, 100, 500, 1000)),
    geom_point(data = full_df, 
               aes(x=obtiene_agua_redpublica_tuberia_adentro_use, y = ALRI5_5_plus_100kpop, 
                   group=period, color=period), alpha=0.5)
  )
)

agua_tuberia_cf <- ggscatterstats(
  data = full_df,
  x = obtiene_agua_redpublica_tuberia_adentro_use,
  y = cf,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Gets water via public system through pipes indoors (%)", # label for x axis
  ylab = "Clean fuel (%)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Piped water", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    ggplot2::scale_y_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    geom_point(data = full_df, 
               aes(x=obtiene_agua_redpublica_tuberia_adentro_use, y=cf, 
                   group=period, color=period), alpha=0.4)
  )
)

agua_tuberia_alri
agua_tuberia_cf

obtiene_agua_alri_mod <- glm(ALRI5_5_plus ~ offset(log(under5population)) + 
                               obtiene_agua_redpublica_tuberia_adentro_use + period + canton_id_universal, 
                                  full_df,  family=quasipoisson())

obtiene_agua_cf_mod <- lm(cf ~ obtiene_agua_redpublica_tuberia_adentro_use + period + canton_id_universal, full_df)

broom::tidy(obtiene_agua_alri_mod)[2,]
broom::tidy(obtiene_agua_cf_mod)[2,]


# Hygiene index ------

hygiene_index_alri <- ggscatterstats(
  data = full_df,
  x = hygiene_index,
  y = ALRI5_5_plus_100kpop,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Hygiene index", # label for x axis
  ylab = "under-5 LRI mortality\n(per 100,000 under-5 population)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Hygiene index", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(),
    scale_y_log10(breaks = c(0.1, 1, 10, 50, 100, 500, 1000)),
    geom_point(data = full_df, 
               aes(x=hygiene_index, y = ALRI5_5_plus_100kpop, 
                   group=period, color=period), alpha=0.5)
  )
)

hygiene_index_cf <- ggscatterstats(
  data = full_df,
  x = hygiene_index,
  y = cf,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Hygiene index", # label for x axis
  ylab = "Clean fuel (%)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Hygiene index", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(),
    ggplot2::scale_y_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    geom_point(data = full_df, 
               aes(x=hygiene_index, y=cf, 
                   group=period, color=period), alpha=0.4)
  )
)


hygiene_index_alri_mod <- glm(ALRI5_5_plus ~ offset(log(under5population)) + 
                                hygiene_index + period + canton_id_universal, 
                             full_df,  family=quasipoisson())

hygiene_index_cf_mod <- lm(cf ~ hygiene_index + period + canton_id_universal, full_df)

broom::tidy(hygiene_index_alri_mod)[2,]
broom::tidy(hygiene_index_cf_mod)[2,]




# electricidad ------

electricidad_alri <- ggscatterstats(
  data = full_df,
  x = electricidad_no,
  y = ALRI5_5_plus_100kpop,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "No grid electricity (%)", # label for x axis
  ylab = "under-5 LRI mortality\n(per 100,000 under-5 population)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "No grid electricity", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    scale_y_log10(breaks = c(0.1, 1, 10, 50, 100, 500, 1000)),
    geom_point(data = full_df, 
               aes(x=electricidad_no, y = ALRI5_5_plus_100kpop, 
                   group=period, color=period), alpha=0.5)
  )
)

electricidad_cf <- ggscatterstats(
  data = full_df,
  x = electricidad_no,
  y = cf,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "No grid electricity (%)", # label for x axis
  ylab = "Clean fuel (%)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "No grid electricity", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    ggplot2::scale_y_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    geom_point(data = full_df, 
               aes(x=electricidad_no, y=cf, 
                   group=period, color=period), alpha=0.4)
  )
)

electricidad_alri
electricidad_cf

electricidad_alri_mod <- glm(ALRI5_5_plus ~ offset(log(under5population)) + 
                                electricidad_no + period + canton_id_universal, 
                              full_df,  family=quasipoisson())

electricidad_cf_mod <- lm(cf ~ electricidad_no + period + canton_id_universal, full_df)

broom::tidy(electricidad_alri_mod)[2,]
broom::tidy(electricidad_cf_mod)[2,]




# Indigenous language spoken: at all ------

idioma_indigena_any_alri <- ggscatterstats(
  data = full_df,
  x = Idioma_indigena_any,
  y = ALRI5_5_plus_100kpop,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Indigenous language spoken at home (%)", # label for x axis
  ylab = "under-5 LRI mortality\n(per 100,000 under-5 population)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Indigenous language spoken at home", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    scale_y_log10(breaks = c(0.1, 1, 10, 50, 100, 500, 1000)),
    geom_point(data = full_df, 
               aes(x=Idioma_indigena_any, y = ALRI5_5_plus_100kpop, 
                   group=period, color=period), alpha=0.5)
  )
)

idioma_indigena_any_cf <- ggscatterstats(
  data = full_df,
  x = Idioma_indigena_any,
  y = cf,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Indigenous language spoken at home (%)", # label for x axis
  ylab = "Clean fuel (%)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Indigenous language spoken at home", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    ggplot2::scale_y_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    geom_point(data = full_df, 
               aes(x=Idioma_indigena_any, y=cf, 
                   group=period, color=period), alpha=0.4)
  )
)

# idioma_indigena2_alri
# idioma_indigena2_cf


# Indigenous language spoken: self ------

idioma_indigena_self_alri <- ggscatterstats(
  data = full_df,
  x = Idioma_indigena_self,
  y = ALRI5_5_plus_100kpop,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Indigenous language spoken at home (%)", # label for x axis
  ylab = "under-5 LRI mortality\n(per 100,000 under-5 population)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Indigenous language spoken at home", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    scale_y_log10(breaks = c(0.1, 1, 10, 50, 100, 500, 1000)),
    geom_point(data = full_df, 
               aes(x=Idioma_indigena_self, y = ALRI5_5_plus_100kpop, 
                   group=period, color=period), alpha=0.5)
  )
)

idioma_indigena_self_cf <- ggscatterstats(
  data = full_df,
  x = Idioma_indigena_self,
  y = cf,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Indigenous language spoken at home (%)", # label for x axis
  ylab = "Clean fuel (%)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Indigenous language spoken at home", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    ggplot2::scale_y_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    geom_point(data = full_df, 
               aes(x=Idioma_indigena_self, y=cf, 
                   group=period, color=period), alpha=0.4)
  )
)


idioma_alri_mod <- glm(ALRI5_5_plus ~ offset(log(under5population)) + 
                               Idioma_indigena_self + period + canton_id_universal, 
                             full_df,  family=quasipoisson())

idioma_cf_mod <- lm(cf ~ Idioma_indigena_self + period + canton_id_universal, full_df)

broom::tidy(idioma_alri_mod)[2,]
broom::tidy(idioma_cf_mod)[2,]



# Indigenous language spoken: parents ------

idioma_indigena_parents_alri <- ggscatterstats(
  data = full_df,
  x = Idioma_indigena_parents,
  y = ALRI5_5_plus_100kpop,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Indigenous language spoken at home (%)", # label for x axis
  ylab = "under-5 LRI mortality\n(per 100,000 under-5 population)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Indigenous language spoken at home", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    scale_y_log10(breaks = c(0.1, 1, 10, 50, 100, 500, 1000)),
    geom_point(data = full_df, 
               aes(x=Idioma_indigena_parents, y = ALRI5_5_plus_100kpop, 
                   group=period, color=period), alpha=0.5)
  )
)

idioma_indigena_parents_cf <- ggscatterstats(
  data = full_df,
  x = Idioma_indigena_parents,
  y = cf,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Indigenous language spoken at home (%)", # label for x axis
  ylab = "Clean fuel (%)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Indigenous language spoken at home", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    ggplot2::scale_y_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    geom_point(data = full_df, 
               aes(x=Idioma_indigena_parents, y=cf, 
                   group=period, color=period), alpha=0.4)
  )
)




# Wall material: hormigon ------

pared_alri <- ggscatterstats(
  data = full_df,
  x = material_pared_hormigon_bloque_ladrillo,
  y = ALRI5_5_plus_100kpop,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "High-quality wall material (%)", # label for x axis
  ylab = "under-5 LRI mortality\n(per 100,000 under-5 population)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "High-quality wall material", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    scale_y_log10(breaks = c(0.1, 1, 10, 50, 100, 500, 1000)),
    geom_point(data = full_df, 
               aes(x=material_pared_hormigon_bloque_ladrillo, y = ALRI5_5_plus_100kpop, 
                   group=period, color=period), alpha=0.5)
  )
)

pared_cf <- ggscatterstats(
  data = full_df,
  x = material_pared_hormigon_bloque_ladrillo,
  y = cf,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "High-quality wall material (%)", # label for x axis
  ylab = "Clean fuel (%)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "High-quality wall material", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    ggplot2::scale_y_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    geom_point(data = full_df, 
               aes(x=material_pared_hormigon_bloque_ladrillo, y=cf, 
                   group=period, color=period), alpha=0.4)
  )
)

pared_alri
pared_cf

pared_alri_mod <- glm(ALRI5_5_plus ~ offset(log(under5population)) + 
                        material_pared_hormigon_bloque_ladrillo + period + canton_id_universal, 
                             full_df,  family=quasipoisson())

pared_cf_mod <- lm(cf ~ material_pared_hormigon_bloque_ladrillo + period + canton_id_universal, full_df)

broom::tidy(pared_alri_mod)[2,]
broom::tidy(pared_cf_mod)[2,]



# Floor material: tierra ------

piso_alri <- ggscatterstats(
  data = full_df,
  x = material_piso_tierra,
  y = ALRI5_5_plus_100kpop,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Floor material is dirt (%)", # label for x axis
  ylab = "under-5 LRI mortality\n(per 100,000 under-5 population)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Floor is dirt", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    scale_y_log10(breaks = c(0.1, 1, 10, 50, 100, 500, 1000)),
    geom_point(data = full_df, 
               aes(x=material_piso_tierra, y = ALRI5_5_plus_100kpop, 
                   group=period, color=period), alpha=0.5)
  )
)

piso_cf <- ggscatterstats(
  data = full_df,
  x = material_piso_tierra,
  y = cf,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Floor material is dirt (%)", # label for x axis
  ylab = "Clean fuel (%)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Floor material is dirt", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    ggplot2::scale_y_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    geom_point(data = full_df, 
               aes(x=material_piso_tierra, y=cf, 
                   group=period, color=period), alpha=0.4)
  )
)

piso_alri
piso_cf

piso_alri_mod <- glm(ALRI5_5_plus ~ offset(log(under5population)) + 
                        material_piso_tierra + period + canton_id_universal, 
                      full_df,  family=quasipoisson())

piso_cf_mod <- lm(cf ~ material_piso_tierra + period + canton_id_universal, full_df)

broom::tidy(piso_alri_mod)[2,]
broom::tidy(piso_cf_mod)[2,]


# Floor material: nicer ------

piso_nicer_alri <- ggscatterstats(
  data = full_df,
  x = material_piso_nicer,
  y = ALRI5_5_plus_100kpop,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Floor material is nicer (%)", # label for x axis
  ylab = "under-5 LRI mortality\n(per 100,000 under-5 population)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Floor is nicer", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    scale_y_log10(breaks = c(0.1, 1, 10, 50, 100, 500, 1000)),
    geom_point(data = full_df, 
               aes(x=material_piso_nicer, y = ALRI5_5_plus_100kpop, 
                   group=period, color=period), alpha=0.5)
  )
)

piso_nicer_cf <- ggscatterstats(
  data = full_df,
  x = material_piso_nicer,
  y = cf,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Floor material is nicer (%)", # label for x axis
  ylab = "Clean fuel (%)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Floor material is nicer", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    ggplot2::scale_y_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    geom_point(data = full_df, 
               aes(x=material_piso_nicer, y=cf, 
                   group=period, color=period), alpha=0.4)
  )
)


piso_nicer_alri_mod <- glm(ALRI5_5_plus ~ offset(log(under5population)) + 
                       material_piso_nicer + period + canton_id_universal, 
                     full_df,  family=quasipoisson())

piso_nicer_cf_mod <- lm(cf ~ material_piso_nicer + period + canton_id_universal, full_df)

broom::tidy(piso_nicer_alri_mod)[2,]
broom::tidy(piso_nicer_cf_mod)[2,]



# Roof material: high quality ------

techo_hormigon_alri <- ggscatterstats(
  data = full_df,
  x = material_techo_hormigon_losa_cemento_use,
  y = ALRI5_5_plus_100kpop,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Roof material is high quality (%)", # label for x axis
  ylab = "under-5 LRI mortality\n(per 100,000 under-5 population)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Roof material: high quality", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    scale_y_log10(breaks = c(0.1, 1, 10, 50, 100, 500, 1000)),
    geom_point(data = full_df, 
               aes(x=material_techo_hormigon_losa_cemento_use, y = ALRI5_5_plus_100kpop, 
                   group=period, color=period), alpha=0.5)
  )
)

techo_hormigon_cf <- ggscatterstats(
  data = full_df,
  x = material_techo_hormigon_losa_cemento_use,
  y = cf,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Roof material is high quality (%)", # label for x axis
  ylab = "Clean fuel (%)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Roof material: high quality", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    ggplot2::scale_y_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    geom_point(data = full_df, 
               aes(x=material_techo_hormigon_losa_cemento_use, y=cf, 
                   group=period, color=period), alpha=0.4)
  )
)

techo_alri_mod <- glm(ALRI5_5_plus ~ offset(log(under5population)) + 
                       material_techo_hormigon_losa_cemento_use + period + canton_id_universal, 
                     full_df,  family=quasipoisson())

techo_cf_mod <- lm(cf ~ material_techo_hormigon_losa_cemento_use + period + canton_id_universal, full_df)

broom::tidy(techo_alri_mod)[2,]
broom::tidy(techo_cf_mod)[2,]


# Roof material: asbestos ------

techo_asbestos_alri <- ggscatterstats(
  data = full_df,
  x = material_techo_asbesto_use,
  y = ALRI5_5_plus_100kpop,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Roof material is asbestos (%)", # label for x axis
  ylab = "under-5 LRI mortality\n(per 100,000 under-5 population)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Roof material: asbestos", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    scale_y_log10(breaks = c(0.1, 1, 10, 50, 100, 500, 1000)),
    geom_point(data = full_df, 
               aes(x=material_techo_asbesto_use, y = ALRI5_5_plus_100kpop, 
                   group=period, color=period), alpha=0.5)
  )
)

techo_asbestos_cf <- ggscatterstats(
  data = full_df,
  x = material_techo_asbesto_use,
  y = cf,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Roof material is asbestos (%)", # label for x axis
  ylab = "Clean fuel (%)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Roof material: asbestos", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    ggplot2::scale_y_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    geom_point(data = full_df, 
               aes(x=material_techo_asbesto, y=cf, 
                   group=period, color=period), alpha=0.4)
  )
)

techo_a_alri_mod <- glm(ALRI5_5_plus ~ offset(log(under5population)) + 
                          material_techo_asbesto_use + period + canton_id_universal, 
                        full_df,  family=quasipoisson())

techo_a_cf_mod <- lm(cf ~ material_techo_asbesto_use + period + canton_id_universal, full_df)

broom::tidy(techo_a_alri_mod)[2,]
broom::tidy(techo_a_cf_mod)[2,]

full_df$material_techo_notnice <- 1 - (full_df$material_techo_asbesto_use + full_df$material_techo_hormigon_losa_cemento_use)
full_df$material_techo_nice <- full_df$material_techo_asbesto_use + full_df$material_techo_hormigon_losa_cemento_use

techo_n_alri_mod <- glm(ALRI5_5_plus ~ offset(log(under5population)) + 
                          material_techo_nice + period + canton_id_universal, 
                      full_df,  family=quasipoisson())

techo_n_cf_mod <- lm(cf ~ material_techo_nice + period + canton_id_universal, full_df)

broom::tidy(techo_n_alri_mod)[2,]
broom::tidy(techo_n_cf_mod)[2,]



# Roof material: higher quality ------

techo_higher_alri <- ggscatterstats(
  data = full_df,
  x = material_techo_nicer_all,
  y = ALRI5_5_plus_100kpop,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Roof material is higher quality (%)", # label for x axis
  ylab = "under-5 LRI mortality\n(per 100,000 under-5 population)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Roof material: higher quality materials", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    scale_y_log10(breaks = c(0.1, 1, 10, 50, 100, 500, 1000)),
    geom_point(data = full_df, 
               aes(x=material_techo_nicer_all, y = ALRI5_5_plus_100kpop, 
                   group=period, color=period), alpha=0.5)
  )
)

techo_higher_cf <- ggscatterstats(
  data = full_df,
  x = material_techo_nicer_all,
  y = cf,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Roof materials are higher quality (%)", # label for x axis
  ylab = "Clean fuel (%)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Roof material: higher quality materials", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    ggplot2::scale_y_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    geom_point(data = full_df, 
               aes(x=material_techo_nicer_all, y=cf, 
                   group=period, color=period), alpha=0.4)
  )
)

techo_nicer_alri_mod <- glm(ALRI5_5_plus ~ offset(log(under5population)) + 
                              material_techo_nicer_all + period + canton_id_universal, 
                            full_df,  family=quasipoisson())

techo_nicer_cf_mod <- lm(cf ~ material_techo_nicer_all + period + canton_id_universal, full_df)

broom::tidy(techo_nicer_alri_mod)[2,]
broom::tidy(techo_nicer_cf_mod)[2,]



# Household materials index ------

materials_index_alri <- ggscatterstats(
  data = full_df,
  x = materials_index,
  y = ALRI5_5_plus_100kpop,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Household building materials index", # label for x axis
  ylab = "under-5 LRI mortality\n(per 100,000 under-5 population)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Household building materials index", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    scale_x_continuous(),
    scale_y_log10(breaks = c(0.1, 1, 10, 50, 100, 500, 1000)),
    geom_point(data = full_df, 
               aes(x=materials_index, y=ALRI5_5_plus_100kpop, 
                   group=period, color=period), alpha=0.5)
  )
)

materials_index_cf <- ggscatterstats(
  data = full_df,
  x = materials_index,
  y = cf,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Househoold building materials index)", # label for x axis
  ylab = "Clean fuel (%)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Househould building materials index", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    scale_x_continuous(),
    ggplot2::scale_y_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    geom_point(data = full_df, 
               aes(x=materials_index, y=cf, 
                   group=period, color=period), alpha=0.4)
  )
)


materials_index_alri_mod <- glm(ALRI5_5_plus ~ offset(log(under5population)) + 
                          materials_index + period + canton_id_universal, 
                        full_df,  family=quasipoisson())

materials_index_cf_mod <- lm(cf ~ materials_index + period + canton_id_universal, full_df)

broom::tidy(materials_index_alri_mod)[2,]
broom::tidy(materials_index_cf_mod)[2,]



# Vaccine: BCG ------

bcg_alri <- ggscatterstats(
  data = full_df,
  x = BCG_use,
  y = ALRI5_5_plus_100kpop,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "BCG vaccine (%)", # label for x axis
  ylab = "under-5 LRI mortality\n(per 100,000 under-5 population)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "BCG vaccine", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    scale_y_log10(breaks = c(0.1, 1, 10, 50, 100, 500, 1000)),
    geom_point(data = full_df, 
               aes(x=BCG_use, y = ALRI5_5_plus_100kpop, 
                   group=period, color=period), alpha=0.5)
  )
)

bcg_cf <- ggscatterstats(
  data = full_df,
  x = BCG_use,
  y = cf,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "BCG vaccine (%)", # label for x axis
  ylab = "Clean fuel (%)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "BCG vaccine", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    ggplot2::scale_y_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    geom_point(data = full_df, 
               aes(x=BCG_use, y=cf, 
                   group=period, color=period), alpha=0.4)
  )
)

bcg_alri
bcg_cf


bcg_alri_mod <- glm(ALRI5_5_plus ~ offset(log(under5population)) + 
                          BCG_use + period + canton_id_universal, 
                        full_df,  family=quasipoisson())

bcg_cf_mod <- lm(cf ~ BCG_use + period + canton_id_universal, full_df)

broom::tidy(bcg_alri_mod)[2,]
broom::tidy(bcg_cf_mod)[2,]


# Vaccine: DPT3  ------

DPT3_carnet_alri <- ggscatterstats(
  data = full_df,
  x = DPT3_use,
  y = ALRI5_5_plus_100kpop,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Diptheria, pertusis, tetanus vaccines (%)", # label for x axis
  ylab = "under-5 LRI mortality\n(per 100,000 under-5 population)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Diptheria, pertusis, tetanus vaccines", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    scale_y_log10(breaks = c(0.1, 1, 10, 50, 100, 500, 1000)),
    geom_point(data = full_df, 
               aes(x=DPT3_use, y = ALRI5_5_plus_100kpop, 
                   group=period, color=period), alpha=0.5)
  )
)

DPT3_carnet_cf <- ggscatterstats(
  data = full_df,
  x = DPT3_use,
  y = cf,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Diptheria, pertusis, tetanus vaccines (%)", # label for x axis
  ylab = "Clean fuel (%)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Diptheria, pertusis, tetanus vaccines", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    ggplot2::scale_y_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    geom_point(data = full_df, 
               aes(x=DPT3_use, y=cf, 
                   group=period, color=period), alpha=0.4)
  )
)

DPT3_carnet_alri
DPT3_carnet_cf

dpt3_alri_mod <- glm(ALRI5_5_plus ~ offset(log(under5population)) + 
                       DPT3_use + period + canton_id_universal, 
                    full_df,  family=quasipoisson())

dpt3_cf_mod <- lm(cf ~ DPT3_use + period + canton_id_universal, full_df)

broom::tidy(dpt3_alri_mod)[2,]
broom::tidy(dpt3_cf_mod)[2,]



# Vaccine: SAR carnet ------

SAR_carnet_alri <- ggscatterstats(
  data = full_df,
  x = SAR_use,
  y = ALRI5_5_plus_100kpop,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Measles vaccine (%)", # label for x axis
  ylab = "under-5 LRI mortality\n(per 100,000 under-5 population)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Measles vaccine", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    scale_y_log10(breaks = c(0.1, 1, 10, 50, 100, 500, 1000)),
    geom_point(data = full_df, 
               aes(x=SAR_use, y = ALRI5_5_plus_100kpop, 
                   group=period, color=period), alpha=0.5)
  )
)

SAR_carnet_cf <- ggscatterstats(
  data = full_df,
  x = SAR_use,
  y = cf,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Measles vaccine (%)", # label for x axis
  ylab = "Clean fuel (%)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Measles vaccine", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    ggplot2::scale_y_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    geom_point(data = full_df, 
               aes(x=SAR_use, y=cf, 
                   group=period, color=period), alpha=0.4)
  )
)


sar_alri_mod <- glm(ALRI5_5_plus ~ offset(log(under5population)) + 
                      SAR_use + period + canton_id_universal, 
                     full_df,  family=quasipoisson())

sar_cf_mod <- lm(cf ~ SAR_use + period + canton_id_universal, full_df)

broom::tidy(sar_alri_mod)[2,]
broom::tidy(sar_cf_mod)[2,]



# Vaccine: POLIO carnet ------

Polio_carnet_alri <- ggscatterstats(
  data = full_df,
  x = Polio_use,
  y = ALRI5_5_plus_100kpop,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Polio vaccine (%)", # label for x axis
  ylab = "under-5 LRI mortality\n(per 100,000 under-5 population)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Polio vaccine", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    scale_y_log10(breaks = c(0.1, 1, 10, 50, 100, 500, 1000)),
    geom_point(data = full_df, 
               aes(x=Polio_use, y = ALRI5_5_plus_100kpop, 
                   group=period, color=period), alpha=0.5)
  )
)

Polio_carnet_cf <- ggscatterstats(
  data = full_df,
  x = Polio_use,
  y = cf,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Polio vaccine (%)", # label for x axis
  ylab = "Clean fuel (%)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Polio vaccine", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    ggplot2::scale_y_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    geom_point(data = full_df, 
               aes(x=Polio_use, y=cf, 
                   group=period, color=period), alpha=0.4)
  )
)

polio_alri_mod <- glm(ALRI5_5_plus ~ offset(log(under5population)) + 
                      Polio_use + period + canton_id_universal, 
                    full_df,  family=quasipoisson())

polio_cf_mod <- lm(cf ~ Polio_use + period + canton_id_universal, full_df)

broom::tidy(polio_alri_mod)[2,]
broom::tidy(polio_cf_mod)[2,]

# Vaccine: pcv3 carnet ------

pcv3_carnet_alri <- ggscatterstats(
  data = full_df,
  x = pcv3_use,
  y = ALRI5_5_plus_100kpop,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "PCV-3 vaccine (%)", # label for x axis
  ylab = "under-5 LRI mortality\n(per 100,000 under-5 population)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Pneumococcal conjugate vaccine - 3", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    scale_y_log10(breaks = c(0.1, 1, 10, 50, 100, 500, 1000)),
    geom_point(data = full_df, 
               aes(x=pcv3_use, y = ALRI5_5_plus_100kpop, 
                   group=period, color=period), alpha=0.5)
  )
)

pcv3_carnet_cf <- ggscatterstats(
  data = full_df,
  x = pcv3_use,
  y = cf,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "PCV-3 vaccine (%)", # label for x axis
  ylab = "Clean fuel (%)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Pneumococcal conjugate vaccine - 3", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    ggplot2::scale_y_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    geom_point(data = full_df, 
               aes(x=pcv3_use, y=cf, 
                   group=period, color=period), alpha=0.4)
  )
)

pcv3_alri_mod <- glm(ALRI5_5_plus ~ offset(log(under5population)) + 
                        pcv3_use + 
                       period +
                       canton_id_universal
                     , 
                      full_df,  family=quasipoisson())

pcv3_cf_mod <- lm(cf ~ pcv3_use + 
                    period +
                    canton_id_universal, 
                  full_df)

broom::tidy(pcv3_alri_mod)[2,]
broom::tidy(pcv3_cf_mod)[2,]



# Vaccine: vaccine index ------

vaccine_index_carnet_alri <- ggscatterstats(
  data = full_df,
  x = vaccine_index,
  y = ALRI5_5_plus_100kpop,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Vaccine index", # label for x axis
  ylab = "under-5 LRI mortality\n(per 100,000 under-5 population)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Vaccine index", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    # ggplot2::scale_x_continuous(labels=scales::percent_format(), 
    #                             breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
    #                             limits = c(0, 1)),
    scale_y_log10(breaks = c(0.1, 1, 10, 50, 100, 500, 1000)),
    geom_point(data = full_df, 
               aes(x=vaccine_index, y = ALRI5_5_plus_100kpop, 
                   group=period, color=period), alpha=0.5)
  )
)

vaccine_index_carnet_cf <- ggscatterstats(
  data = full_df,
  x = vaccine_index,
  y = cf,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Vaccine index", # label for x axis
  ylab = "Clean fuel (%)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Vaccine index", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    # ggplot2::scale_x_continuous(labels=scales::percent_format(), 
    #                             breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
    #                             limits = c(0, 1)),
    ggplot2::scale_y_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    geom_point(data = full_df, 
               aes(x=vaccine_index, y=cf, 
                   group=period, color=period), alpha=0.4)
  )
)

vaccine_index_alri_mod <- glm(ALRI5_5_plus ~ offset(log(under5population)) + 
                       vaccine_index + period + canton_id_universal, 
                     full_df,  family=quasipoisson())

vaccine_index_cf_mod <- lm(cf ~ vaccine_index + period + canton_id_universal, full_df)

broom::tidy(vaccine_index_alri_mod)[2,]
broom::tidy(vaccine_index_cf_mod)[2,]

# average maternal age ------

mage_alri <- ggscatterstats(
  data = full_df_x,
  x = mage_mean,
  y = ALRI5_5_plus_100kpop,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Average age (years)", # label for x axis
  ylab = "under-5 LRI mortality\n(per 100,000 under-5 population)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Average maternal age at birth", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    # ggplot2::scale_x_continuous(labels=scales::percent_format(), 
    #                             breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
    #                             limits = c(0, 1)),
    scale_y_log10(breaks = c(0.1, 1, 10, 50, 100, 500, 1000)),
    geom_point(data = full_df_x, 
               aes(x=mage_mean, y = ALRI5_5_plus_100kpop, 
                   group=period, color=period), alpha=0.5)
  )
)

mage_cf <- ggscatterstats(
  data = full_df_x,
  x = mage_mean,
  y = cf,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Maternal age (years)", # label for x axis
  ylab = "Clean fuel (%)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Average maternal age", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    # ggplot2::scale_x_continuous(labels=scales::percent_format(),
    #                             breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
    #                             limits = c(0, 1)),
    ggplot2::scale_y_continuous(labels=scales::percent_format(),
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    geom_point(data = full_df_x, 
               aes(x=mage_mean, y=cf, 
                   group=period, color=period), alpha=0.4)
  )
)


mage_alri_mod <- glm(ALRI5_5_plus ~ offset(log(under5population)) + 
                      mage_mean + period + canton_id_universal, 
                    full_df_x,  family=quasipoisson())

mage_cf_mod <- lm(cf ~ mage_mean + period + canton_id_universal, full_df_x)

broom::tidy(mage_alri_mod)[2,]
broom::tidy(mage_cf_mod)[2,]


# any ANC visits ------

ANC_alri <- ggscatterstats(
  data = full_df,
  x = ANC_use,
  y = ALRI5_5_plus_100kpop,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Any ANC visits (%)", # label for x axis
  ylab = "under-5 LRI mortality\n(per 100,000 under-5 population)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "ANC visits", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    scale_y_log10(breaks = c(0.1, 1, 10, 50, 100, 500, 1000)),
    geom_point(data = full_df, 
               aes(x=ANC_use, y = ALRI5_5_plus_100kpop, 
                   group=period, color=period), alpha=0.5)
  )
)

ANC_cf <- ggscatterstats(
  data = full_df,
  x = ANC_use,
  y = cf,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Any ANC visits (%)", # label for x axis
  ylab = "Clean fuel (%)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Any ANC visits", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    ggplot2::scale_x_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    ggplot2::scale_y_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    geom_point(data = full_df, 
               aes(x=ANC_use, y=cf, 
                   group=period, color=period), alpha=0.4)
  )
)

ANC_alri
ANC_cf

ANC_alri_mod <- glm(ALRI5_5_plus ~ offset(log(under5population)) + 
                      ANC_use + period + canton_id_universal, 
                              full_df,  family=quasipoisson())

ANC_cf_mod <- lm(cf ~ ANC_use + period + canton_id_universal, full_df)

broom::tidy(ANC_alri_mod)[2,]
broom::tidy(ANC_cf_mod)[2,]


# ANC visits median ------

ANC_visits_alri <- ggscatterstats(
  data = full_df,
  x = ANC_visits_median_use,
  y = ALRI5_5_plus_100kpop,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Median number of ANC visits", # label for x axis
  ylab = "under-5 LRI mortality\n(per 100,000 under-5 population)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Median number of ANC visits", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    # ggplot2::scale_x_continuous(labels=scales::percent_format(), 
    #                             breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
    #                             limits = c(0, 1)),
    scale_y_log10(breaks = c(0.1, 1, 10, 50, 100, 500, 1000)),
    geom_point(data = full_df, 
               aes(x=ANC_visits_median_use, y = ALRI5_5_plus_100kpop, 
                   group=period, color=period), alpha=0.5)
  )
)

ANC_visits_cf <- ggscatterstats(
  data = full_df,
  x = ANC_visits_median_use,
  y = cf,
  # grouping.var = period,
  type = "robust", # type of test that needs to be run
  conf.level = 0.95, # confidence level
  xlab = "Median number of ANC visits", # label for x axis
  ylab = "Clean fuel (%)", # label for y axis
  line.color = "black", # changing regression line color line
  title = "Median ANC visits", # title text for the plot
  
  ggtheme = theme_bw(), # choosing a different theme
  ggstatsplot.layer = FALSE, # turn off ggstatsplot theme layer
  marginal.type = "density", # type of marginal distribution to be displayed
  xfill = "#42B540B2", # color fill for x-axis marginal distribution
  yfill = "#00468BB2", # color fill for y-axis marginal distribution
  xalpha = 0.6, # transparency for x-axis marginal distribution
  yalpha = 0.6, # transparency for y-axis marginal distribution
  point.alpha = 0.2, 
  messages = FALSE, # turn off messages and notes
  ggplot.component = list(
    # ggplot2::scale_x_continuous(labels=scales::percent_format(), 
    #                             breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
    #                             limits = c(0, 1)),
    ggplot2::scale_y_continuous(labels=scales::percent_format(), 
                                breaks = c(0, 0.20, 0.40, 0.60, 0.80, 1),
                                limits = c(0, 1)),
    geom_point(data = full_df, 
               aes(x=ANC_visits_median_use, y=cf, 
                   group=period, color=period), alpha=0.4)
  )
)


ANC_visits_alri_mod <- glm(ALRI5_5_plus ~ offset(log(under5population)) + 
                             ANC_visits_median_use + period + canton_id_universal, 
                    full_df,  family=quasipoisson())

ANC_visits_cf_mod <- lm(cf ~ ANC_visits_median_use + period + canton_id_universal, full_df)

broom::tidy(ANC_visits_alri_mod)[2,]
broom::tidy(ANC_visits_cf_mod)[2,]


# 4. Correlation between covariates ---------

covs_corr <- ggstatsplot::ggcorrmat(
  
  data = full_df_x %>% 
    mutate(vaccine_index = -vaccine_index) %>% 
    dplyr::select(cf, percent_rural_corrected_ratio,
                    electricidad_no,
                  material_techo_nicer_all,
                    # material_techo_hormigon_losa_cemento_use,
                    material_pared_hormigon_bloque_ladrillo,
                  material_piso_nicer,
                  # material_piso_tierra,
                    materials_index,
                    
                    obtiene_agua_redpublica_tuberia_adentro_use,
                    # servicio_higenico_redpublica_use,
                  # elimina_servidas_red_pozo,
                  elimina_servidas_red_pozo_letrina,
                    servicio_ducha_excluso,
                    elimina_basura_service,
                  hygiene_index,
                    
                    # Literate, 
                  Literate_women,
                    # InSchool, 
                  InSchool_girls,
                    # Idioma_indigena_any, 
                  Idioma_indigena_self, 
                  # Idioma_indigena_parents,
                    
                    BCG_use,
                    DPT3_use,
                    SAR_use,
                    Polio_use,
                    pcv3_use,
                    vaccine_index,
                    
                    mage_mean, 
                  ANC_use,
                    ANC_visits_median_use), 
  cor.vars.names = c("Clean Fuel Use",
                     "Rural",
                     "Not Grid Electrified",
                     
                     "Roof Material: Nicer",
                     "Wall Material: Nicest",
                     "Floor Material: Nicer",
                     "Materials Index",
                     
                     "Municipal Water Piped Inside",
                     "Modern Toilet with Waste Removal,\nCesspool/Septic Tank,\nLatrine",
                     "Private Shower",
                     "Municipal Trash Removal",
                     "Household Hygiene Index",
                     
                     "Adult Female Literacy",
                     "Girls under 18 School Attendance",
                     "Indigenous Language Spoken at Home",
                     "u-5 Tuberculosis Vaccine",
                     "u-5 Diptheria, Pertussis,\nTetanus Vaccine",
                     "u-5 Measles Vaccine",
                     "u-5 Polio Vaccine",
                     "u-5 Pneumonia Conjugate (3) Vaccine",
                     "u-5 Vaccine Index",
                     "Average age of mothers at delivery",
                     "Antenatal Care",
                     "Median Antenatal Care Visits"),
  sig.level = 1,
  title = "Correlations between all potential confounding variables",
  lab.size = 3.5,
  caption = ""
  
) 

